* (c) M. Ciccarelli. Rats version 7.0
*
** program to do forecast using previously estimated components.
** compare the performance with RW AR PHIL. Reproduce Table 8.
**
** the model is the usual stock-watson
*
** y_h(t+h) = c + a(L)*y(t) + d(L)*F(t)+ e(t)
*
** where notation is as in the paper.
** In particular:
*
** 1. y_h(t+h) = (400/h)log(P(t)/P(t-h))
*
** 2. y(t) = 400*log(P(t)/P(t-1))
*
*  change sample at line 168 for Phillips curve results
*
** The model is compared withthree benchmarks
*
** a. AR
** b. Phillips curve
** c. RW
*
*

OPEN DATA forecast_DATA.xls    ;* SA data
CALENDAR 1960 1 4
ALL 2008:02
DATA(FORMAT=XLS,ORG=COLUMNS) 1960:01 2008:02
close data
comp lags1 = 1
comp seas = 1


** deseasonalise data
if seas ==1
{
X11(GRADUATE) CPIQU2 / CPIQU2
X11(GRADUATE) CPIQAT / CPIQAT
X11(GRADUATE) CPIQLU / CPIQLU
X11(GRADUATE) CPIQIT / CPIQIT
X11(GRADUATE) CPIQPT / CPIQPT
X11(GRADUATE) CPIQNL / CPIQNL
X11(GRADUATE) CPIQSE / CPIQSE
X11(GRADUATE) CPIQBE / CPIQBE
X11(GRADUATE) CPIQDE / CPIQDE
X11(GRADUATE) CPIQDK / CPIQDK
X11(GRADUATE) CPIQES / CPIQES
X11(GRADUATE) CPIQFI / CPIQFI
X11(GRADUATE) CPIQuk / CPIQUK
X11(GRADUATE) CPIQFR / CPIQFR
X11(GRADUATE) CPIQGR / CPIQGR
X11(GRADUATE) CPIQIE / CPIQIE
X11(GRADUATE) CPIQUS / CPIQUS
X11(GRADUATE) CPIQJP / CPIQJP
X11(GRADUATE) CPIQAU / CPIQAU
X11(GRADUATE) CPIQCA / CPIQCA
X11(GRADUATE) CPIQNZ / CPIQNZ
X11(GRADUATE) CPIQNO / CPIQNO
X11(GRADUATE) CPIQCH / CPIQCH
}
end if

comp n = 23
comp nsteps = 8
comp lags = 0

dec vec[series] fore_diff(n)
dec vec[series] fore_infl(n)
dec vec[series] national(n)
dec vec[series] infla(n)
dec vec[series] infla1(n)
dec vec[series] infla_steps(n)
dec vec[series] diff_infla(n)
dec vec[series] idio(n)

dec vect[series] ip(n) m3(n)

dec vec[series] diffe(n)
dec vec[series] diffe_sq(n)
dec vec[series] diffe_abs(n)

dec vec[series] diffe_naiv(n)
dec vec[series] diffe_naiv_sq(n)
dec vec[series] diffe_naiv_abs(n)

dec vec rmse_far(n) mad_far(n)
dec vec rmse_naiv(n) mad_naiv(n)
dec vec rmse_ar1(n) mad_ar1(n)
dec vec rmse_far1(n) mad_far1(n)

dec vec[series] diffe_ar1(n)
dec vec[series] diffe_ar1_sq(n)
dec vec[series] diffe_ar1_abs(n)

dec vec[series] diffe_far1(n)
dec vec[series] diffe_far1_sq(n)
dec vec[series] diffe_far1_abs(n)

dec vec[series] diffe_phil(n)
dec vec[series] diffe_phil_sq(n)
dec vec[series] diffe_phil_abs(n)

dec vec rmse_phil(n) mad_phil(n)
dec vec[series] fore_diff_phil(n)
dec vec[series] fore_infl_phil(n)

dec vec[series] fore_diff_ar1(n)
dec vec[series] fore_infl_ar1(n)

dec vec[series] fore_diff_rw1(n)
dec vec[series] fore_infl_rw1(n)

dec vec[series] fore_diff_far1(n)
dec vec[series] fore_infl_far1(n)
dec vec[series] resid(n)
dec rec[rec] ratios(3,nsteps)

do i = 1,3
do j = 1,nsteps
dim ratios(i,j)(n,2)
end
end

comp lag_bic = 2;*  max lag to search with bic criterion
comp lagopt = 1;* 0 = no optimal choice of lag
comp lag_fact1 = 4
comp lag_fact2 = 2
comp lag_fix = 1
COMP AVER = 0

dec vec[integer] opt_lag(n)
dec vec bic(lag_bic)

source(noecho) princomp.src
source(noecho) uforeerr.src

dofor steps = 1 4 8

disp steps
*
do i = 1,n
set infla(i) = 400.*log(([series] i)/(([series] i){1}))
set infla_steps(i) = (400./steps)*log(([series] i)/(([series] i){steps}))
end

do i = 1,n
set national(i) = infla_steps(i)-infla(i){steps}
end

comp J=1;
dofor i = M3QU2 M3QUS	M3QCA	M3QUK	M3QJP	M3QDE	M3QFR	M3QIT	M3QAT	M3QBE	M3QDK	M3QFI	M3QGR	M3QIE	M3QDE M3QPT	M3QES	M3QSE	M3QNL	M3QAU	M3QNZ	M3QNO	M3QCH
set m3(j) = log(i{0}/i{1})*400
comp J=J+1;
end dofor

comp J=1;
dofor i = IPQU2	IPQUS	IPQCA	IPQUK	IPQJP	IPQDE	IPQFR	IPQIT	IPQAT	IPQBE	IPQDK	IPQFI	IPQGR	IPQIE	IPQLU	IPQPT	IPQES	IPQSE	IPQNL	IPQAU	IPQNZ	IPQNO	IPQCH
set ip(j) = log(i{0}/i{1})*400.0
comp J=J+1;
end dofor

set dcp = log(comprice/comprice{1})*400

comp beg = 1960:3
comp tempend = 1995:4
comp end = 2008:2    ;* FOR THE PHILLIPS CURVE CHANGE TO 2007:2
comp fin = end-tempend

smpl beg end

clear forint FORE_INFL FORE_INFL_RW1 FORE_INFL_AR1 FORE_INFL_PHIL

comp inum = 0

 DO ti = steps,fin
*
  comp count = 1
   dofor jj = infla(1) infla(2) infla(3) infla(4) infla(5) infla(6) infla(7) infla(8) $
    infla(9) infla(10) infla(11) infla(12) infla(13) infla(14) infla(15) infla(16)$
    infla(17) infla(18) infla(19) infla(20) infla(21) infla(22) infla(23)
    stat(noprint) jj beg tempend+ti-steps
    set infla1(count) = (([series] jj)-%mean)/sqrt(%variance)
    comp count = count+1
   end

  @PRINCOMP(Ncomp=2) beg tempend+ti-steps pcomp
  # infla1(1) infla1(2) infla1(3) infla1(4) infla1(5) infla1(6) infla1(7) infla1(8) $
    infla1(9) infla1(10) infla1(11) infla1(12) infla1(13) infla1(14) infla1(15) infla1(16)$
    infla1(17) infla1(18) infla1(19) infla1(20) infla1(21) infla1(22) infla1(23)

IF AVER==0
{
   set factor1 beg tempend+ti-steps = pcomp(1)
}
ELSE
{
  set average beg tempend+ti-steps = 0.0

   do i = 1,n
   set average beg tempend+ti-steps = average+(1./n)*infla(i)
   end
   set factor1 beg tempend+ti-steps = average
}

do i = 1,n

  if lagopt==1
  {
  do maxlag=1,lag_bic
   linreg(noprint) INFLA_steps(i) beg tempend+ti-steps
   # constant infla(i){steps to steps+maxlag-1} factor1{steps to steps+lag_fact1-1}
   compute bic(maxlag) = log(%rss/%nobs)+log(%nobs)*%nreg/%nobs
  end do

  do maxl = 1,lag_bic
   comp max = %minvalue(bic)
   if max == bic(maxl)
   comp opt_lag(i) = fix(maxl)
  end
  }
  else
  {
   comp opt_lag(i) = lag_fix
  }

  equation eq infla_steps(i)
  # constant infla(i){steps to steps+opt_lag(i)-1} factor1{steps to steps+lag_fact1-1}
  linreg(equation=eq,noprint) INFLA_steps(i) beg tempend+ti-steps
  forecast 1 steps tempend+ti-steps+1
  # eq forint
  set fore_infl(i) tempend+ti tempend+ti = forint
  clear forint

  end do i

 comp inum = inum+1

 END do TI

*** COMPUTE STATISTICS to compare with competitors
*

 do i = 1,n
   @uForeErrors(noprint) infla_steps(i) fore_infl(i) tempend+steps end
   comp rmse_far(i) = %%FERRMSE
 end

do i = 1,n
set diffe(i) tempend+steps end = fore_infl(i)-infla_steps(i)
set diffe_sq(i) tempend+steps end = (fore_infl(i)-infla_steps(i))**2.
set diffe_abs(i) tempend+steps end = abs(fore_infl(i)-infla_steps(i))
end

************************
** compute statistics NAIVE BENCHMARK

smpl beg end

clear forint

comp inum = 0

DO ti = steps,fin

do i = 1,n
linreg(define=eq,noprint) INFLA_steps(i) beg tempend+ti-steps
# infla(i){steps}
RESTRICT(create,noprint,define=eq)  1 resids1
# 1
# 1.0 1
forecast 1 steps tempend+ti-steps+1
# eq forint
set fore_infl_rw1(i) tempend+ti tempend+ti = forint
clear forint
end
comp inum = inum+1

END do TI

 do i = 1,n
   @uForeErrors(noprint) infla_steps(i) fore_infl_rw1(i) tempend+steps end
   comp rmse_naiv(i) = %%FERRMSE
 end

do i = 1,n
set diffe_naiv(i) tempend+steps end = fore_infl_rw1(i)-infla_steps(i)
set diffe_naiv_sq(i) tempend+steps end = (fore_infl_rw1(i)-infla_steps(i))**2.
set diffe_naiv_abs(i) tempend+steps end = abs(fore_infl_rw1(i)-infla_steps(i))
end

**************

** compute statistics SIMPLE AR(p) FORECAST


clear forint
comp inum = 0

DO ti = steps,fin

do i = 1,n

  if lagopt==1
  {
  do maxlag=1,lag_bic
   linreg(noprint) INFLA_steps(i) beg tempend+ti-steps
   # constant infla(i){steps to steps+maxlag-1}
   compute bic(maxlag) = log(%rss/%nobs)+log(%nobs)*%nreg/%nobs

  end do

  do maxl = 1,lag_bic
   comp max = %minvalue(bic)
   if max == bic(maxl)
   comp opt_lag(i) = fix(maxl)
  end
  }
  else
  {
   comp opt_lag(i) = lag_fix
  }
equation eq infla_steps(i)
# constant infla(i){steps to steps+opt_lag(i)-1}
linreg(equation=eq,noprint) infla_steps(i) beg tempend+ti-steps
forecast 1 steps tempend+ti-steps+1
# eq forint
set fore_infl_ar1(i) tempend+ti tempend+ti = forint
clear forint
end

comp inum = inum+1
END do TI

 do i = 1,n
   @uForeErrors(noprint) infla_steps(i) fore_infl_ar1(i) tempend+steps end
   comp rmse_ar1(i) = %%FERRMSE
 end

do i = 1,n
set diffe_ar1(i) tempend+steps end = fore_infl_ar1(i)-infla_steps(i)
set diffe_ar1_sq(i) tempend+steps end = (fore_infl_ar1(i)-infla_steps(i))**2.
set diffe_ar1_abs(i) tempend+steps end = abs(fore_infl_ar1(i)-infla_steps(i))
end

**************

** compute statistics SIMPLE Phillips curve FORECAST
*
clear forint

comp inum = 0

DO ti = steps,fin
*disp ti

do i = 1,n

  if lagopt==1
  {
  do maxlag=1,lag_bic
   linreg(noprint) INFLA_steps(i) beg tempend+ti-steps
   # constant infla(i){steps to steps+maxlag-1} dcp{steps to steps+maxlag-1} ip(i){steps to steps+maxlag-1} m3(i){steps to steps+maxlag-1}
   compute bic(maxlag) = log(%rss/%nobs)+log(%nobs)*%nreg/%nobs

  end do

  do maxl = 1,lag_bic
   comp max = %minvalue(bic)
   if max == bic(maxl)
   comp opt_lag(i) = fix(maxl)
  end
  }
  else
  {
   comp opt_lag(i) = lag_fix
  }
equation eq infla_steps(i)
# constant infla(i){steps to steps+opt_lag(i)-1} dcp{steps to steps+opt_lag(i)-1} ip(i){steps to steps+opt_lag(i)-1} m3(i){steps to steps+opt_lag(i)-1}
linreg(equation=eq,noprint) infla_steps(i) beg tempend+ti-steps
forecast 1 steps tempend+ti-steps+1
# eq forint
set fore_infl_phil(i) tempend+ti tempend+ti = forint
clear forint
end
comp inum = inum+1

END do TI

 do i = 1,n
   @uForeErrors(noprint) infla_steps(i) fore_infl_phil(i) tempend+steps end
   comp rmse_phil(i) = %%FERRMSE
 end

do i = 1,n
set diffe_phil(i) tempend+steps end = fore_infl_phil(i)-infla_steps(i)
set diffe_phil_sq(i) tempend+steps end = (fore_infl_phil(i)-infla_steps(i))**2.
set diffe_phil_abs(i) tempend+steps end = abs(fore_infl_phil(i)-infla_steps(i))
end

**************

*** TEST FOR EQUAL PREDICTIVE ACCURACY
** compare with RANDOM WALK FORECAST

do i = 1,n
compute ratios(1,steps)(i,1) = rmse_far(i)/rmse_naiv(i)
*set F_CW = diffe_naiv_sq(i) - (diffe_sq(i)- ((fore_infl_rw1(i) - fore_infl(i))**2.))
set F_CW = diffe_naiv_sq(i) - diffe_sq(i)
linreg(noprint,robust,lag=steps-1,damp=1) F_CW tempend+steps end
# constant
if %tstats(1).gt.1.282
{
comp ratios(1,steps)(i,2) = 1.0
}
else
{
comp ratios(1,steps)(i,2) = 0.0
}
clear F_CW
end

** compare with SIMPLE AR(p) FORECAST

do i = 1,n
compute ratios(2,steps)(i,1) = rmse_far(i)/rmse_ar1(i)
set F_CW = diffe_ar1_sq(i) - diffe_sq(i)
*set F_CW =  diffe_ar1_sq(i) - (diffe_sq(i)- ((fore_infl_ar1(i) - fore_infl(i))**2.))
linreg(noprint,robust,lag=steps-1,damp=1) F_CW tempend+steps end
# constant
if %tstats(1).gt.1.282
{
comp ratios(2,steps)(i,2) = 1.0
}
else
{
comp ratios(2,steps)(i,2) = 0.0
}
clear F_CW
end

** compare with Augmented Phillips curve FORECAST

do i = 1,n
compute ratios(3,steps)(i,1) = rmse_far(i)/rmse_phil(i)
set F_CW = diffe_phil_sq(i)-diffe_sq(i)
linreg(noprint,robust,lag=steps-1,damp=1) F_CW tempend+steps end
# constant
if %tstats(1).gt.1.282
{
comp ratios(3,steps)(i,2) = 1.0
}
else
{
comp ratios(3,steps)(i,2) = 0.0
}
clear F_CW
end

DO I = 1,N
clear diffe(i) diffe_sq(i) diffe_abs(i)
clear diffe_naiv(i) diffe_naiv_sq(i) diffe_naiv_abs(i)
clear diffe_phil(i) diffe_phil_sq(i)
clear diffe_ar1(i) diffe_ar1_sq(i)
END

end dofor


** write here the ratio of the RMSE as in table 5 at horizons=4 8 steps.

if end == 2007:2
{
disp 'results for Phillips curve columns in table 5'
disp '4-step ahead'
wri ratios(3,4) ;* PHil
disp '8-step ahead'
wri ratios(3,8)
}
if end == 2008:2
{
disp 'results for RW columns in table 5'
disp '4-step ahead'
wri ratios(1,4) ;* RW
disp '8-step ahead'
wri ratios(1,8)

disp 'results for AR columns in table 5'
disp '4-step ahead'
wri ratios(2,4) ;* AR
disp '8-step ahead'
wri ratios(2,8)
}
end if





